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ABSTRACT 

We present Smoothed Particle Hydrodynamics simulations of protoclusters including 
the effects of the stellar winds from massive stars. Using a particle-injection method, 
we investigate the effect of structure in close proximity to the wind sources and the 
short-timescale influence of winds on protoclusters. We find that structures such as 
disks and gaseous filaments have a strong collimating effect on winds. By a different 
technique of injecting momentum from point sources into our simulations, we compare 
the large-scale and long-term effects of isotropic and intrinsically-collimated winds 
on protoclusters and find them to be similar, although the collimated winds take 
longer to exert a significant influence. We find that both types of wind are able 
to dramatically slow the global star formation process, but that the timescale on 
which they can expel significant quantities of mass from the cluster is rather long 
(approaching ten freefall times). Clusters may then experience rapid star formation 
very early in their lifetimes, before switching to a mode where gas is gradually ex- 
pelled, while star formation proceeds much more slowly over many freefall times. This 
complicates any conclusions regarding slow star formation derived from measuring 
the star-formation efficiency per freefall time. We find that estimates of the efficacy 
of winds in dispersing clusters derived simply from the total wind momentum flux 
may not be very reliable. This is due to material being expelled from deep within 
stellar potential wells, ofen to velocities well in excess of the cluster escape velocity, 
and also to the loss of momentum flux through holes in the gas distribtuion. Winds 
have little effect on the accretion-driven stellar initial mass function (IMF) except 
at the very high-mass end, where they serve to prevent some of the most massive 
objects accreting more material. Feedback does not result in the formation of further 
massive stars through the monolithic collapse of massive cores. We also find that 
the morphology of the gas, the rapid motions of the wind sources and the action of 
large-scale accretion flows prevent the formation of bubble-like structures. These 
effects may make it difficult to discern the influence of winds on very young clusters. 
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1 INTRODUCTION 

Elucidating the effects of stellar feedback on embedded star 
clusters is a problem of key importance in astrophysics. 
~ 90% of newborn stellar systems do not survive their 
embed ded phase to become open clusters (|Lada fc Ladal 
(2003)). In addition, star formation at the scale of molecular 
cloud s and embedded c l usters is not a very efficient process 
(e.g. iFukui fc Miaunol (|l99ll )). There are two possible 



* E-mail: jim@ig.cas.cz (JED) 



reasons for these observations. If star-forming complexes 
are initially gravitationally bound (as is usually assumed), 
some process must expel most of the molecular gas before it 
can be incorporated into stars, while also usually disrupting 
the system. Alternatively, if such complexes are not initially 
bound, the low star formation efficiency and low survival 
probability are a natural consequence of this, although 
bound stellar systems c an be left behin d in s ome cases 
l|Clark fc Bonnelll (|2004T ). IClark fc Bormell <|2005l )h 

If star-forming systems are initially bound, an obvious 
candidate for a means of unbinding them is feedback from 



2 J. E. Dale, LA. Bonnell 



their own stars. In embedded clusters massive enough to 
form O-type stars (i.e. those with masses > 10 3 Mq, the 
ionizing radiation, winds and eventual supernova (SN) ex- 
plosions of these stars are likely to dominate the combined 
effect of the (much more numerous) lower mass objects. 
Photoionization and winds are particularly important in 
this respect, since they act for the duration of the O-star's 
lifetime, whereas SNe cannot occur for at least 5 Myr 
after the onset of star formation. This corresponds to 
the freefall time of a molecular cloud of number density 
~ 10 2 cm' 3 , approximately the minimum density required 
to for m molecular gas able to cool efficiently and form 
stars (|Hartmann et all {2001)). If it is the case that star 
formation in, and disruption of, molecular clo uds occurs 
on tim escales of around the local freefa ll time (Elmcgrccn 
|200d ). but see iKrumholz fc Tanl (|2007T )), SNe may occur 
too late, particularly as the density of a molecular cloud 
generally increases considerably (with a corresponding 
decrease in the freefall time) before any massive stars are 
able to form. 

In this paper, we consider the effects of the winds from 
massive stars on the cluster in which they form. In reality, 
winds and photoionization act in concert and it is not clear 
whi ch will dominate, nor if one wil l help or hinder the other 
(see ICapriotti fc Kozminskil (|2001) for a discussion of this 
subject). In this paper, we model winds alone - we leave a 
comparative study including winds and ionization to a later 
date. 

The interaction of a stellar wind with the inter- 
stellar medium (ISM) has been modeled analytically 
and numerically by several authors, althoug h usua lly in 
one-dimensio n al ge ometries. ISteigman et al.l (1 19751 ) and 
IWeaver et alj l| 19771 ) introduced the idea of dividing the 
wind bubble into several concentric zones. The innermost 
zone in their models consists of freely-flowing wind and the 
outermost zone of undisturbed ISM. The number and con- 
tents of the zones in be t ween is the subject of considerable 
debate. ISteigman et al.l (|l975l ) have a single zone between 
the free-flowing wind and the undisturbed ISM, consisting 
of a shocked shell driven by the momentum of the wind. It is 
not difficult to show that this assumption produces a wind 
bubb le whose radius R(i) increases with time as R ( t) oc t 2 



(e.g. lOstriker fc McKed ll98Sl ; lLamers fc Cassinellil Il999h ) 



IWeaver et all (|l977n introduce two zones in this region - an 
inner shell of hot shocked wind and an outer shell of shocked 
ISM separated by a contact discontinuity. The difference 
between the models lies in Steigman et al's assumption 
that the shocked wind can cool and that the bubble is 
momentum-driven, as opposed to Weaver et al's contention 
that the shocked wind cannot cool and that the bubble 
is thus pressure-driven (resulting in a radial evolution 
of the form R(t) oc i§). ICapriotti fc Kozminskil (|200ll ) 
discuss at some length the validity of these assumptions, 
which essentially rests on whether the contact discontinuity 
between the shocked wind and the shocked ISM is stable, 
which governs how likely mixing and subsequ e nt co oling 
of these phases is. Analytical work b y Vishniad (ll983T l and 
simulations by iGarcia-Segura et al.l ( 19961 ) suggest that 
the discontinuity may not be stable and the shocked wind 
may therefore be able to cool efficiently. In this case, the 



wind bubble is driven purely by ram pressure and R(t) oc 1 2 . 



2 NUMERICAL METHODS 

The code used in these calculation s is a hybrid SPH/N- 
body code described in iBate et"aL (1995). Hydrodynamic 
equations are solved using the SPH formalism and stars are 
modeled as collisionless point masses which are able to ac- 
crete gas particles. We have included the effects of stellar 
winds emanating from point masses in two ways. In the first 
instance, we model the wind directly by injecting low-mass 
SPH particles. The second method models the sink particles 
as momentum sources and distributes the momentum flux 
amongst the SPH gas particles near the wind sources. 



2.1 Particle winds 

The first method of handling stellar winds involves injecting 
low-mass SPH particles from the chosen sink-particles. In 
order to establish a continuous wind from the ejection of 
discrete particles, a series of shells are emitted at regular 
time intervals At. Each shell consists of n w ind (typically 
300-500) SPH particles which are chosen randomly from 
6141 positions equally spaced over the surface of a sphere, 
with random rotations between sucessive shells to ensure 
that particles do not cluster along any particular directions. 
The shells radial width is determined from the wind 
velocity, v^, and the time interval between each shell's 
ejection: Ar = VooAt. The r-coordinates of the particles 
injected into the innermost shell are chosen randomly in 
the range [rmin , rmin + Ar] , where r m in must be outside 
the accretion radius of the sink particle. The mass of the 
wind particles are then set to provide the desired mass loss 
rate. Provided that there are a reasonable number (> 10) 
shells between the wind source and the nearest regular gas 
particle, this procedure produces an intrinsically smooth 
wind. However, owing to to the very high wind terminal 
velocity, this dictates that the timestep At on which winds 
must be updated be very small, which makes the code 
comparatively slow. 

In order that the interaction of wind particles with 
regular gas particles at the contact surface is well modelled, 
the smoothing lengths of the wind particles should be 
similar to those of the regular gas particles, so that both 
types of particle should have 50 — 100 neighbours of either 
particle type. Once a contact surface with an inner skin of 
wind particles is established, this is achieved reasonably 
well by the SPH code without assistance. However, when 
the wind is first launched, to avoid particle intepenetration 
at the contact surface, the smoothing lengths of the wind 
particles must be set so that they are similar to those of 
the regular gas particles with which they will first come 
into contact. We first define a smoothing length h ncal and 
a distance i? nC ar representative of the regular particles 
nearest the wind source from the initial gas distibution. 
Since the smoothing lengths of the wind particles evolve 



in the spherical wind as h oc 



„2/3 



the wind-particles' 



smoothing lengths can be set so that /i w ind(i?ncar) ~ linear) 
by an appropriate combination of are then set initially from 
a combination of n w i n( j,Ar (determined by the choice of 
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(a) Shock radius, uniform medium 
(p =constant). Solid line is the analytic 
fit R(t) oc r 1 / 2 , circles arc the results 
from the simulation. 



(b) Shock ^-displacement, slab with p oc 
J/ - 1 . Solid line is the analytic fit Y(t) oc 
f 2 / 3 , crosses are results from the simula- 
tion. 




(c) Particle plot of particles within 2h of 
the x-plane (y on the vertical axis, z on 
the horizontal axis), slab with p oc \y\~ 1 . 
Blue crosses are regular gas particles, red 
particles are wind particles. 



Figure 1. Results of tests of the particle-injection wind code. Quantities are in code units. 



At) and r m i n . We set r m i n to be just outside the accretion 
radius of the wind source (10 _4 pc in all the simulations in 
this paper) and choose n w i n d to be 300-500 to ensure that 
individual wind shells are reasonably smooth, and use this 
to constrain At. In the simulations presented here, these 
quantities are estimated at the beginning of each run and 
kept constant. It would improve the code to set their values 
dynamically, but we find in the simple test cases below 
and the more challenging simulations presented later that 
particle interpenetration is minimal and that the shock 
surfaces are well modelled. Once a shock has been well 
established, particle interpenetration is less of a problem, 
since the the dense particles at the contact surface prevent 
it. Another improvement to the method, not implemented 
here, would require that the smoothing lengths are calcu- 
lated only considering particles of the same type (i.e. wind 
particles must have ~ 50 wind particle neighbours, and 
regular particles must have ~ 50 regular neighbours). When 
particles of very different properties (masses, temperatures 
etc) are present in the same simulation, the evolution of the 
smoothing lengths should consider the number of particles 
with similar properties. In the simulations involving highly 
inhomogenous gas distributions, this would help limit 
potential problems due to lower density regions being 
under-resolved and thus the artificial creation of voids by 
the winds. 

The particle-injection wind code was tested in two 
simple cases: (i) a single wind source in a uniform medium 
(ii) a single source embedded in a slab whose density profile 
follows p oc The analytic solution for the radius of the 

wind b ubble in the former c ase, R (t) oc t 1 ^ 2 , can be found 
in, e.g, lOstriker fc McKeel (|l988). In the latter case, we 
consider a disc perpendicular to the y-axis and of constant 
size plouging through the density distribution and sweeping 
up mass proportional to J Q p(y)dy. The flux absorbed 
by the disc is proportiona to y~ 2 , so the size of the wind 
bubble in the y-direction as Y(t) oc t 2 ^ 3 . In Figures |l(a)| 
and 1(b) we show the results of these tests compared to the 
appropriate analytic functions. Figure 1(c) shows a particle 



plot from the simulation of the wind source embedded 
in a slab, and demonstrates that there is no interparticle 
penetration. We conclude that the particle-injection code 
performs well in these simple test cases. 

In Figure [2j we show the results of the evolution of a 
wind bubble in a highly inhomogenous medium. We use 
the standard SPH viscosity but with the values a = 2 and 
(5 = 4 in order to ensure the shock is adequately captured. 
The wind particles are hot (10 4 K, corresponding to an 
ionised wind) and are assumed to behave isothermally. 
We again find that there is no interparticle penetration. 
The great difficulty in using SPH particles to model the 
stellar winds is that their dynamical timescales are much 
shorter than that of the other gas particles present in the 
simulation. This means that the evolution can be followed 
only for relatively short timescales and long term effects 
cannot be investigated. This is the principal reason why the 
second method, described below, has also been developed. 



2.2 Momentum— injection winds 

The second method is to inject momentum into the simu- 
lation and distribute it amongst the regular SPH gas par- 
ticles. This has the advantage that fast-moving wind par- 
ticles, with corresponding small timesteps, need not be in- 
troduced into the simulation. This allowed us to follow the 
evolution of the wind bubbes for much longer times, so that 
we could study the effects of winds on the early evolution of 
protoclusters. 

In line with the work of IVishniad (Il983l ) and 
iGarcia-Segura et al.l (|l996l ) , we consider winds driven purely 
by ram pressure (although it must be admitted that this 
choice is also partly motivated by computational necessity, 
since simulation of pressure-driven winds requires detailed 
radiative transfer calculations which are beyond the scope 
of this work). 

Each source is first assigned a momentum flux. The SPH 
particles around the source over which this flux is to be dis- 
tributed are then identified. We define around each wind 
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Figure 2. The outflowing gas from a massive star's wind is shown 
(panels 2-4) as it is collimated by the filamentary structure in the 
surrounding cluster gas (far left panel). Images are 0.05 pc 10 4 
au) on a side. The column density of the gas is shown from 1.5 
to 400 g cm -2 in the top left panel while the column density of 
the wind extends from 10 -5 to 0.02 g cm -2 . 



source a 'working face' (by analogy with a coal face) con- 
sisting of the minimum number of SPH particles which are 
required to shadow or shield all the particles in the simula- 
tion from the wind source. In order to shield a given particle, 
a working-face particle must cut a line joining that particle's 
centre to the wind source. A momentum packet leaving the 
source must then either strike a working-face particle (be- 
fore striking any other particles) or strike no particles at all 
and exit the cloud In which case the working face does not 
cover all of 4ir steradians). We illustrate this idea in Figure 

El 

Once the working face is defined, the momentum flux 
received by each working face particle must be determined. 
Since most working face particles are partially obscured by 
others, the fraction of the source's momentum flux absorbed 
by each one is not simply dictated by the solid angle the par- 
ticle subtends at the source. What is required instead is the 
unobscured solid angle subtended by each working face par- 
ticle. In principle, this could be calculated geometrically by 
considering the particles as overlapping circles (seen from 
the point of view of the souce). However, this problem be- 
comes very complicated when there are many such mutually 
overlapping circles of different sizes, as will be the case in a 
complex density distribution with SPH particles of different 
sizes at different distances from the source. We therefore do 
not explicitly calculate the unobscured solid angle presented 
to the source by the working face particles. Instead, we es- 
timate the flux they each receive using a method analagous 
to that employed by Monte Carlo radiative transfer codes, 
although the problem at hand is much simpler. We imagine 
that the wind momentum flux is divided up into a number 
of packets which are emitted in random directions from the 
wind source. Each packet either strikes an SPH particle or 



does not - it is likely in the highly-inhomogeneous gas distri- 
butions existing in molecular clouds that some stellar wind 
exits through a hole in the cloud and strikes the interstellar 
medium outside, but we do not concern ourselves with this 
process. We assume that if a momentum packet strikes an 
SPH particle, all of the momentum in the packet is absorbed 
by that particle. If a packet does not strike any particles, we 
ignore it (although the momentum of the packet is still deb- 
ited from the source's total momentum flux). 

Up to ten passes are performed in which 100 x n/ ace 
packets are launched, where n/ ace is the number of particles 
in the working face. The wind code exits when either all 
particles in the working face have been hit at least once, or 
ten passes are complete. Note that, owing to partial shad- 
owing by other working face particles, the proportion of hits 
any member of the working face is likely to receive is highly 
variable, with a few particles very close to the wind source 
likely to be struck many times, while those far from the 
source which are partially shielded by other particles may 
only receive a few impacts. It is this issue that necissitates 
the use of the Monte Carlo method, since the momentum 
flux recieved by a working face particle is not simply. The 
momentum flux received by each particle in the working 
face is then determined by multiplying the wind flux from 
the source by the fraction of emitted packets hitting that 
particle. 

We considered three test cases: (i) a single wind source 
in a uniform medium (ii) a single source in a spherically 
symmetric density distribution following p oc r^ 1 (iii) a sin- 
gle source embedded in a slab w hose density profile follow s 
p cx \y\ 1/2 . Referring again to lOstriker fc McKeei l|l988l) . 
we find that the wind bubbles in the first two cases fol- 
low R(t) oc i 1 / 2 and R(t) oc t 2/l3 . In the third case, we again 
consider a disc of constant size moving in the {/-direction 
and obtain the size of the wind bubble in the {/-direction 
as Y(t) oc i 4 / 7 . We show the results of these tests in Figure 
[4] and compare them to appropriate analytic solutions. In 
each case, the extent of the wind bubble was followed by 
tracking the positions of the densest gas particles, i.e. those 
in the shock driven by the wind. The slight jaggedness in 
the plot of the simulated data in Figure 4(c) is a result of 



the necessity of using only a few tens of particles to accu- 
rately track the tips of the wind bubble in the positive and 
negative {/-directions. We find that the code reproduces the 
expected behaviour very well. 



3 INTRINSIC AND EXTRINSIC 
COLLIMATION 

Given the complex structures of molecular clouds, it is likely 
that an isotropic wind emitted by a star will be collimated 
to some degree by density inhomogeneities near the star. 
This effect has been observed in simulations o f the emis- 
sion o f photoionizing radiation from O-stars by Dale et al.l 
i|2005l ). and we observe it in our particle-injection wind sim- 
ulations, as detailed in Section 4, where we show that small- 
scale structure very close to the wind sources can collimate 
outflows into bipolar or even unipolar morphologies. For ex- 
ample, in Figure[7]the presence of a circumstellar disc results 
in a wind bubble formed that resembles an hourglass, rather 
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(a) Shock radius, uniform medium 
(p =constant). Dashed line is the analytic 
fit R(t) oc t 1 / 2 , solid line is the result 
from the simulation. 
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(b) Shock radius, medium with p oc r _1 . 
Dashed line is the analytic fit R(t) oc t 2 / 3 , 
solid line is the result from the simulation. 
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(c) Shock ^-displacement, slab with p oc 
\y\~ 1 / 2 . Dashed line is the analytic fit 
Y(t) oc t 4 / 7 , solid line is the result from 
the simulation. 



Figure 4. Results of tests of the momentum— injection wind code. 
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Figure 3. Illustration of the working face. Note that the particle 
marked with a X symbol is not shielded by the working face 
particle marked with a + as that working face particle does not 
cut the line joining the wind source to the centre of the particle 
marked with X . 



than a sphere. We refer to such small-scale collimation as 
'intrinsic' collimation. 

In our long-timescale protocluster simulations, in which 
we use the momentum-injection technique to model winds, 
it is likely that the winds would be intrinsically collimated 
by small-scale structure close to the wind sources which is 
not well resolved. To investigate the effect this might have 
on the evolution of the protocluster, we include the effects of 
intrinsic collimation by modifying the angular dependence 
of the momentum output. We consider primarily collimation 
due to a circumstellar disc. Although we may not be able 
to resolve the disk itself, our code does calculate the angu- 
lar momentum of sink particles, derived from the angular 



momentum of the material from which the sink originally 
formed, plus that from material accreted later. As the pu- 
tative accretion disk is very likely to have an angular mo- 
mentum vector oriented in the same direction as that of the 
young star around which it orbits, we can use the instan- 
taneous angular momentum vector of each sink particle to 
define the axis of the unresolved accretion disk. 

Once an axis has been defined, the simplest way to 
achieve a polar wind is to modulate the initially-isotropic 
distribution of wind momentum using a function propor- 
tional to the cosine of the polar angle, raised to some power 
n which determines the degree of collimation. There is no 
real physics involved in this - it is simply a convenient de- 
vice by which we can explore the general importance of col- 
limation. Modulating the wind momentum flux in this way 
would result in a lower radiated wind power. We perform 
simulations in which we compare the effects of collimated 
and uncollimated winds. In the case of collimated winds, we 
increase the source momentum fluxes so that the actual ra- 
diated wind power of the collimated and uncollimated wind 
sources is the same. For an arbitrary power n, we achieve 
this by multiplying the wind momentum flux by a normal- 
izing factor F given by 



F = 4tt 



cos"9sin9d9 



(1) 



In our simulations in which we employ anisotropic 
winds, we use n = 5, so that the polar angle angle at 
which the winds flux falls to half its maximum value (which 
gives a measure of the wind opening half-angle) is ~ 29°, 
corresponding to a wind weakly collimated by a disk. 
By comparison with a simulation in which the winds are 
not intrinsically collimated, we can study separately the 
effects of small-scale collimation by circumstellar disks, 
and collimation on larger scales by accretion flows or other 
structures in the gas. We term this collimation on larger 
scales 'extrinsic'. 

The Monte Carlo technique is relatively computation- 
ally expensive, although the algorithm can be parallelised 
efficiently. The overhead from the mometum injection wind 
code in calculations involving a single wind source depends 
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quite strongly on how many particles are in the working 
face. In the three tests detailed above, involving single 
sources and run on a single CPU machine, the fractional 
runtime consumed by the wind code ranged from negligible 
in the early stages of simulations to up to ~ 80% towards 
the end. In the calculations detailed below in Section 5, 
run on 4-CPU machines and in which 6 — 12 sources 
are involved, the simulations with winds active took ~ 5 
times longer to run than those without. Nothwithstanding 
this considerable overhead, the momentum-injection code 
is very much faster than the particle-injection code and 
enables calculations to continue for timescales on the order 
of cluster freefall times so that the efects of winds on early 
cluster evolution, for example gas expulsion and cluster 
unbinding, may be studied. 









Jt 



Figure 5. AitofT projections of the cold cluster gas (left) and 
the hot wind (right) are shown divided into regions close to the 
star (< 500 AU, top) and further away (> 500 AU, bottom). We 
see that the cluster gas impedes the wind in certain directions, 
thereby collimating it into the channels where the column density 
of the cluster gas is smallest. 



4 PARTICLE WINDS AND 

INTRINSICALLY— COLLIMATED 
OUTFLOWS 

We apply our particle-injection method to test the intrinsic 
collimation of isotropic winds. In order to do this we 
make use of the massive star s formed in prev i ous n umer- 
ical simulatio n s rep orted in iBonnell fc Bate! (|2002l ) and 
iBonnell et ail (|2003l ). In the former, we model only the 
central part of the cluster containing 327 stars of total mass 
151 Mq embedded in 371 Mq of gas. The wind is from 
a 30 Mq star located at the centre of the system and is 
surrounded by a remnant filament of gas (see first panel of 
figure [2]). The wind is emitted just beyond the sink-radius 
of 10 AU with a mass outflow rate of 1 x 10~ 5 M yr _1 , 
a wind velocity of 1000 km s _1 and an internal sound 
speed of 10 km s _1 . The evolution was followed for a total 
of 85 years, by which time the wind consisted of 90,000 
particles and there were 77,000 regular gas particles. The 
total outflow distance is only 5600 AU, corresponding to an 
average expansion velocity of ~ 300km s _1 . If taken at face 
value, this would imply that the wind particles are slowed 
by interaction with a quantity of surrounding gas equal to 
only approximately four times their own mass and that the 
interaction of the two fluids is not well resolved. We find in 
fact that the wind particles travel almost ballistically until 
they strike the sorrounding gas, rapidly depositing almost 
all their momentum and coming virtually to a halt. There is 
no interparticle penetration and the interaction of the wind 
with the surrounding gas is well resolved. It is the short 
timesteps required to integrate this process that makes it 
impractical to evolve these calculations for long timescales. 
We also find, however, that the wind particles find and fill 
holes in the surrounding gas. Since low density regions are 
not well resolved in SPH, this process may not always be 
treated accurately. 

Figure [2] shows the evolution of the wind as it encoun- 
ters the surrounding medium. The filamentary structure of 
the surrounding gas acts to impede the progress of the wind, 
which eventually breaks out only in one direction producing 
a unipolar outflow that is highly collimated. This collima- 
tion of a spherical wind due simply to the external medium 
occurs even in the absence of any rotationally supported disc 



in the simulation, being caused instead by the filamentary 
gas structure in which the massive star is embedded. AitofT 
projections shown in Figure [5] show that there is no parti- 
cle interpenetration present in the simulation. Instead, the 
wind is halted when it interacts with the cold cluster gas. It 
is only in directions in which the column density of the sur- 
rounding gas is low that the wind is able to clear a channel 
through which it can escape. 

Incorporating such winds into larger scale simulations 
is also possible although at the expense of no longer follow- 
ing the long-term evolution of the full syste m. We take the 
same i nitial conditions as detailed below (from ffBonnell et al.l 
(2003)), and inject particle-based winds. All stars whose 
mass M t exceeds 10 Mq are considered to be wind sources. 
Each such object is first assigned a wind mass loss rate M 
according to the expression 

M = W~ 5 ( J^!—) M yr _1 . (2) 
V30M©/ ° y K 1 

All sources are taken to have a wind terminal velocity Voo 
of 10 3 km s -1 . Wind particles are injected at the sink radii 
(200 AU) of the six stars whose mass exceeds 10 Mq and 
the evolution is followed for « 440 years (Fig. [6} ■ We note 
that even without any material within 200 AU, the flows 
are likely to emerge as bipolar due to the collimation by 
structures present at larger scales in the cloud. 

Figure [7] shows a higher resolution simulation of the 
wind from the most massive star in the above calculation. In 
this case, we resolve structures down to 50 AU from the star 
while the wind is launched at 20 AU. This allows the resolu- 
tion of circumstellar discs in the simulation (see Clark et. al 
2007) that further collimate the stellar winds. This results 
in better collimation of the outflows into classical bipolar 
morphologies. The bright ring of material in the wind shows 
the interaction with the circumstellar disc that is providing 
the collimation of the wind. 

The above results demonstrate that any spherical wind 
is likely to be collimated by the external environment and of- 
ten by structures on small scales which are difficult to resolve 
in simulations. Unfortunately, using direct particle winds is 
computationally very expensive, so such simulations cannot 
be followed for large fractions of the lifetime of a typical pro- 
tocluster. In addition, the filling of holes in the surrounding 
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Figure 6. The winds from six massive stars (right) are shown as they interact with the cluster environment (left). Initial conditions 
taken from the self-consistent cluster formation simulation of Bonnell et. al (2003). Note that only four winds are visible as two of the 
six sources are in the same binary system. The winds are shown 440 years after they were launched and the box is 0.5 pc on a side. 




Figure 7. The winds from the most massive star from a higher 
resolution simulation where circumstellar structures are resolved 
down to 50 AU. The discs present at these scales are able to colli- 
mate the outflows into bipolar morphologies. The figure shows the 
hot gas present in the wind in a projection where the circumstellar 
disc is nearly edge-on. The image shows a region approximately 
4000 AU on a side. 



gas by the wind particles raises issues of resolution, since 
these regions are by definition not well resolved and may 
not be true holes. These issues are exacerbated by the use 
of two different fluids with different temperatures, particle 
masses and smoothing lengths. One improvement, which we 
do not incorporate here, is to set smoothing lengths based 
on the requirement that a particle's smoothing kernel must 
contain the canonical number of fifty neighbours of particles 
of the same fluid type (i.e. 'wind' or 'regular'). We do not 
attempt this extension here, since the more basic problem of 
the short timesteps required by the particle injection scheme 
drove us to develop the momentum-driven wind method in 
order to explore the larger-scale and longer-term effects of 
feedback from massive stars on the forming protocluster. 



5 WINDS AND PROTOCLUSTER 
EVOLUTION 

The timescales for which we can compute the evolution of 
the particle-driven winds are necessarily short. In order to 
follow the long-term evolution of a forming protocluster, we 
use our momentum-driven winds method to follow the feed- 
back from the most massive stars in a self-consistent simu- 
lation. 



5.1 Initial Conditions 

We take as ou r initi al conditions simulations performed by 
iBonnell et alj |2003l ). These authors studied the evolution 
of a 10 3 Mq molecular cloud, initially 1 pc in diameter with 
a gas temperature of 10 K. The cloud was initially spheri- 
cal, but seeded with a divergenceless Gaussian supersonic 
turbulent velocity field with a power spectrum P(k) oc fc -4 , 
where k is the wavenumber of the velocity perturbation. 
The total kinetic energy of the cloud of the cloud was set 
equal to the modulus of the gravitational potential energy, 
making the cloud marginally bound. The thermal Jeans 
mass of the cloud was 1 Mq and an isothermal equation of 
st at e was used thr o ughou t . 

IBonnell et all (|2003h followed the evolution of the 
protocluster for 4.5 x 10 5 yr, equivalent to approximately 
2.6 freefall times. They observed that star formation in the 
system was hierarchical. They found that the gas initially 
fragmented at several locations, producing a number of 
small high-density clusters after w 2.4 x 10 5 yr (1.4 t//). 
Over the next freefall time or so, these smaller clusters fell 
in towards each other and merged, forming a single large, 
cent rally-condensed clus ter containing m 400 stars. 

IBonnell et al.l (|200&t ) point out, however, that several 
stars whose mass exceeds 10 Mq form by the end of the 
simulation, and their mechanical and radiative feedback 
may be expected to have a dramatic effect on the subse- 
quent evolution of the system. In the calculations presented 
here, we model the mechanical feedback, i.e. the winds from 
these few massive stars. 

We turn on winds at 3.1 x 10 J yr, ~ 1.8 freefall times 
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Figure 8. Column-density map projected along the z-axis of the 
initial conditions for the momentum-injection wind simulations 
presented in this paper. 

into the cluster's evolution. By this stage, the hierarchical 
star formation process has produced five well-defined 
subclusters and there are six stars whose mass exceeds 10 
Mq. The total number of stars is 377. The total stellar 
mass is sa 411 Mq, so that the star formation efficiency at 
this stage of the system's evolution is ps 41%. We show a 
column-density map of the cluster at this epoch in Figure 
[8] Stars are shown as white dots. We perform three runs 
- a control run which is allowed to evolve without the 
action of winds and two 'feedback' runs in which isotropic 
or anisotropic winds are active. We again treat any object 
whose mass exceeds 10 Mq as a wind source. All wind 
sources have mass-loss rates given by Equation [2] and wind 
terminal velocities of 1000 km s -1 . The momentum fluxes 
Mvao and wind luminosities of our sources as a 

function of mass are shown graphically in Figure [5] 

We follow the evolution of the feedback runs and the 
control run for a further 1.3 freefall times until the cluster 
is ~ 3.1 freefall times old. We followed the control run 
for a slightly shorter time, until an age of ~ 2.9 freefall times. 



6 RESULTS II: PROTO-CLUSTER 
EVOLUTION 

In Figure [TU] we plot the evolution with time of the total 
stellar mass and the star-formation efficiency (SFE) in our 
three simulations. We see that the SFE increases approxi- 
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Figure 9. Source wind momentum fluxes and wind luminosities 
as a function of mass as given by Equation [2] 

mately linearly with time in the run without feedback, as 
would be expected for a gravitationally-bound isothermal 
system. By contrast, in both runs without feedback, 
although the SFE also increases linearly with time in the 
early parts of the simulations, the increase begins to tail 
off. The departure from the behaviour of the no-feedback 
calculation is gradual in both cases and occurs later in the 
run with anisotropic winds (at ~ 2.4 t//) than it does in the 
run with spherically-symmetric winds (at ~ 2.2 t//). The 
winds do not induce the formation of many new stars and 
do not produce any large cores which subsequently collapse 
monolithically to form massive stars. Very few new stars 
form in any calculation, so the discrepancy between the 
SFE of the three runs at any given time is almost entirely 
due to different accretion rates of pre-existing stars. 
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Figure 11. Evolution with the time of the masses of the nine most massive objects in the control run (solid lines), compared with the 
evolution of the same objects in the isotropic (dashed line) and anisotropic (dot-dashed line) winds runs. 



This effect is clear when accretion on the the wind 
sources themselves is examined. In Figure 1111 we plot the 
accretion histories of the nine most massive stars in the 
run without feedback and compare their accretion histories 
with the same objects in the wind-affected calculations. In 
all cases, the effect of these stars' winds is to cut the rate at 
which they accrete. Source 3, whose initial mass is ~ 23 Mq, 
in particular accretes almost nothing in the isotropic winds 
run and very little in the anisotropic winds run during the 
1.3 iff (3.3 x 10 5 ) yr duration of our simulation. This is an 
excellent example of the self-regulation of feedback. The 



same effect is seen to varying degrees in all of Sources 1-5. 
Sources 6-9, however, are different in that their masses 
in the winds runs may be larger or smaller than their 
counterparts in the control run. It is also evident from 
Source 7 that accretion onto a given object may be strongly 
increased by the action of nearby winds, since this source 
crosses the 10 Mq threshold much earlier in the anisotropic 
winds run than in the control run. It is clear from Figure [TT] 
that, although the effect of winds on the global accretion 
rate in the cluster may be simply described, the effect on 
the accretion rate of individual objects, whether they are 
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Figure 10. Evolution with time of the total stellar mass and star- 
formation efficiency in the simulation without feedback (solid 
line), the simulation with isotropic wind feedback (dotted line) 
and the simulation with intrinsically-collimated wind feedback 
(dashed line). 



wind sources or not, is much more complicated. 

In Figure 1121 we explore the self-regulation issue by 
plotting the wind luminosity from all stars with M» > 10 
Mq (as computed using Equation [2} as a function of time 
in the wind runs and comparing it to the total flux that 
all such sources in the windless run would have if they 
were turned on at that instant. We see that the wind 
luminosity levels off in the isotropic wind run as accretion 
onto the wind sources slows almost to a halt. In contrast, 
the fictional wind luminosity in the run without feedback 
continues to rise almost linearly as accretion proceeds 
unchecked. The anisotropic winds run lies between these 
extremes, as accretion onto the winds sources is not slowed 
as efficiently. 

Since the accretion rates onto stars are clearly affected 
by winds, we examined the global effect of winds on 
the initial mass functions (IMFs) of our model clusters, 
depicted in Figure [ToH The similarity of the mass functions 
below ~ 8M0 demonstrates that winds have little effect 
on the IMF of low- and intermediate-mass stars in these 
simulations. The IMFs do differ somewhat above ~ 8M0, 
approximately the same range of masses as the objects 
considered to be wind sources (> IOMq). Winds, as also 
shown in Figure 1111 slow or stop the accretion of the 
wind sources themselves and can prevent them evolving 
to higher masses. This is seen in Figure IT3"1 where the two 
most massive sources in the isotropic winds run (the red 
dashed line) lie in the mass range 20 — 3OM0, whereas 
these objects have evolved to be 30 — 40Mq in the run 
without winds (the solid black line). It also appears that a 
few objects in the isotropic winds run have piled up in the 
10 — 2OM0 mass bin, instead of evolving to higher masses 
as in the other runs. The IMF in the isotropic winds run 
thus appears slightly steeper at the high-mass end, but this 



interpretation should be treated with some caution, as it is 
based on ~ 10 objects. The high-mass IMF of the run with 
intrinsically-collimated winds is intermediate between the 
two other runs. As mentioned earlier, winds do not induce 
the formation of may new stars and do not result in the 
formation of massive clumps giving rise to massive stars - 
the massive stars in our calculations acquire almost all of 
their mass through accretion. We find therefore that the 
IMF in these calculations is still generated almost entirely 
by competitive accretion with little interference from winds. 

To explore the effect of the winds on the global 
dynamics of the cluster, we plot in Figure Q3] the evolution 
of the kinetic, gravitational potential and total energies 
of our cluster (the thermal energies in all runs are both 
constant and negligible in comparison to the other forms of 
energy) . The short-timescale fluctuations in the kinetic and 
gravitational energies in all cases are due to close stellar 
interactions. In the run without feedback, we see that there 
is a general increase with time of global kinetic energy, 
but that this is more than counterbalanced by the global 
increase in magnitude of the gravitational energy, so that 
the system becomes more bound as time progresses. This 
behaviour is what one would expect from an initially-bound 
isothermal system. In the two runs in which winds are 
included, we see rather different behaviour. In both cases, 
the kinetic energy also increases with time. Towards the 
end of the run with anisotropic winds, the rate of increase 
becomes greater than in either of the other calculations. 
Additionally, in both simulations including winds, although 
the gravitational energy initially decreases in a similar 
manner to that in the run without feedback, the decrease 
is slowed somewhat. This is due to expulsion of material 
by the stellar winds, resulting in a corresponding drop in 
the gravitational binding energy. The overall result is that 
the total energy in the winds runs levels off and starts to 
become less negative at ~ 3.0t//. Both the isotropic and 
anisotropic winds are unbinding the model cluster. 

In Figure 1151 we plot the evolution with time of the 
quantity of unbound mass (defined as the sum of the masses 
of all sink- or SPH-particles with positive overall energy) 
in our three calculations. 

We see that the quantity of unbound mass is small 
(< 10%) and roughly constant in the run without feedback, 
as might be expected. The run in which isotropic winds 
are acting (the dotted line in Figure shows an unbound 
mass fraction increasing approximately linearly from the 
beginning of the simulation. The increase is not great 
however, with only ~ 15% of the cluster being unbound 
after winds have been acting for over a freefall time. Once 
again, the run including anisotropic winds shows behaviour 
intermediate between the other two calculations. For a 
period of ~ 0.6 freefall times after the ignition of the 
anisotropic wind sources, they have almost no effect on the 
quantity of unbound mass in the cluster, as compared with 
the run in which winds of any kind are absent. However, 
after ~ 0.6 freefall times, the quantity of unbound mass 
in the anisotropic winds run also begins to rise linearly at 
approximately the same rate as in the run with isotropic 
winds, to a similarly modest value of ~ 12% after ~ 1.3 
freefall times. 

A comparison of the morphologies of the densest 
regions of the three model clusters is shown in Figure 1161 
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Surprisingly, the morphologies of the simulations are rather 
similar, so much so that they are difficult to tell apart. 
There is no well-defined wind bubble in either of the two 
runs in which winds are included. This appears to be due 
to the fact that the wind sources are moving around at high 
velocities into and then out of dense gas, rather than sitting 
in one place. They therefore do not have the opportunity to 
build up a single bubble. Accretion into the central regions 
of the cluster also contributes to the inability of the massive 
stars to construct a wind bubble. Instead of a largely radial 
outflow emanating from the most massive stars, gas is being 
driven out of the cluster anisotropically. The only notable 
difference between the run without feedback and the runs 
with winds are that the winds runs have quantities of gas 
at higher densities than exist in the control run. This is due 
to shock compression of some of the gas by the stellar winds. 



Figure 12. Evolution with time of the total instantaneous wind 
luminosities in the control run (solid line, assuming all winds 
sources were suddenly turned on at each instant), in the simula- 
tion with isotropic wind feedback (dotted line) and the simulation 
with intrinsically-collimated wind feedback (dashed line). 
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Figure 13. The IMFs of the control run (solid black line), 
the isotropic winds run (solid red line) and the intrinsically- 
collimated winds run (dotted green line). Masses are given in 
solar masses. 



7 DISCUSSION 

If one thinks of accretion and expulsion from the system 
as two sinks of the remaining gas in our model cluster, one 
can see which of the two is dominant from the ratio of the 
rate at which material is becoming unbound from the sys- 
tem to the star formation rate (SFR). In Figure [T7] we plot 
this ratio as a function of time for our calculations. If the 
ratio is less than one (below the dashed line in the plots), 
gas is being converted to stars faster than it is being ex- 
pelled, in which case feedback is not the dominant factor 
controlling the dynamics of the gas. Conversely, if this ra- 
tio exceeds one, gas is being expelled from the cluster faster 
than it is being converted to stars and feedback is dominant. 
Not surprisingly, in the run without feedback, the ratio is al- 
ways less than one. In the calculations in which we included 
winds, the ratio is less than unity for most of the duration 
of our simulations, although shows a general upward trend 
as the simulations progress (note, however, that neither the 
global accretion rate nor the rate at which material becomes 
unbound are steady functions, so there are considerable de- 
viations from this trend). In both the simulations including 
winds, after winds have been acting for ~ a freefall time, 
the rate at which material is being expelled from the cluster 
generally exceeds the rate at which it is being accreted by 
existing stars or forming new ones, so feedback is in con- 
trol of the gas dynamics. However, from Figures [lOl and [Pol 
we know that this trend is not due to an acceleration in the 
rate at which the system is losing gas, but rather to a dropoff 
in the global gas accretion/star formation rate. ~ 1 freefall 
time after winds are activated (and ~ 3 freefall times into 
the lifetime of the protocluster), in both winds runs, gas is 
being unbound at an approximately constant rate of ~ 50 
Mq per freefall time. Given that, in both cases, there is still 
~ 400 Mq of gas remaining in the cluster in both simula- 
tions and the accretion/star formation rate is tailing off, we 
conclude that feedback is likely to take ~ 8 freefall times to 
clear the system of gas by the action of winds alone. 
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Figure 14. Evolution with time of the global kinetic (solid line), gravitational potential (dot-dashed line) and total (dashed line) energies 
in the runs without feedback (left panel), with isotropic winds (centre panel) and with anisotropic winds (right panel). 
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Figure 15. Evolution with time of the quantity of unbound mass 
in the control run (solid lines), compared with the evolution of 
the same objects in the isotropic (dashed line) and anisotropic 
(dot-dashed line) winds runs. 



8 CONCLUSIONS 

It is clear from these simulations that spherical stellar 
winds can be collimated by structures near the wind 
sources at a wide range of scales. Whether the winds are 
isotropic or intrinsically anisotropic (i.e. collimated by some 
structure very close to the source) has some influence on 
the effectiveness of the winds, but largely to delay the 
point at which winds become important in controlling 
the cluster dynamics and the star-formation process. We 
found that our intrinsically-collimated winds had similar 
effects to the isotropic winds, but delayed by ~ 0.2 freefall 
times, a relatively small fraction of the duration of our 
simulations. In either case, winds can have a strong effect 



on the evolution of protoclusters. 

We have shown that winds can dramatically slow 
the process of conversion of molecular gas into stars. In 
our simulations, this was achieved largely by slowing or 
stopping accretion on existing stars, rather than preventing 
the formation of new stars, since few new stars formed 
in any of our simulations after the point at which wind 
sources were introduced (~ 1.8 freefall times into the 
evolution of the system). In particular accretion onto the 
more massive wind sources themselves, which constitute a 
significant fraction (~ 20%) of the total stellar mass, was 
severely curtailed by isotropic or anisotropic winds. Since 
the accretion rates onto the intermediate- and low-mass 
stars were not so strongly affected, we found that winds 
had little discernible effect on the cluster IMF except at 
the high-mass end, at and above a mass similar to the 
threshold mass we chose to determine whether a star was 
considered to be a wind source or not. It is possible that 
the (slight) changes in the IMF are therefore influenced 
by our choice of this parameter, so we do not attempt to 
draw strong conclusions on this subject, except that winds 
appear to have little influence on the generation of the IMF 
by competitive accretion. Although there do not appear to 
be any instances of star formation induced by winds in our 
calculations, this may only be a consequence of the fact 
that the clusters are only marginally bound and that the 
wind sources are turned on relatively late, after much of 
the bound gas has already fragmented. It s possible that, in 
clusters that are more strongly bound, winds acting earlier 
in the cluster's evolution may encourage fragmentation and 
induce star formation. 

Additionally, the stellar winds from a few stars can 
come to dominate the dynamics of a whole cluster, so that 
the rate at which mass is being expelled from the cluster 
comes to exceed the rate at which mass is being converted 
to stars. In this case, the winds are beginning to unbind the 
cluster and an examination of the evolution of the kinetic, 
gravitational and total energies of our model clusters 
revealed that the action of winds is to stop and even reverse 
the decline in the cluster's total energy. However, the rates 
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at which winds were able to expel gas from the clusters in 
our simulations were relatively low (~ 50 Mq per freefall 
time, or ~ 5 percent of the cluster's total mass per freefall 
time). This amounts to ~ 3 x lO -4 M0yr _1 . As can be seen 
by comparison with Figure [9] this is much less than the 
momentum fluxes of the sources, if the gas is expelled at 
velocities comparable to the cluster escape velocity of a 
few kilometres per second. We identify three reasons for 
this discrepancy. Firstly, gas is generally expelled from 
regions near the massaive stars, where escape velocities are 
- 100km s _1 . Secondly, some of the gas leaving the cluster 
is moving at velocities in excess of the escape velocity. 
Thirdly, due to the inhomogenous gas distribution, some of 
the momentum flux from the sources is simply lost through 
holes in the gas. This result implies that the momentum 
flux of embedded wind sources alone may not be a good 
criteroin for determining whether a given protocluster can 
be dispersed by winds. 

The slow rate of gas expulsion led to an interesting 
situation in which our wind-influenced clusters still pos- 
sessed large quantities of gas which it would take the winds 
nearly ten freefall times to expel, but were converting gas 
to stars very slowly. The time for which the clusters were 
forming stars at a the rate dictated by purely isothermal 
evolution (as observed in the control run without winds) 
would then be only ~ 20 percent of the time for which they 
possessed molecular gas. If generally true, this would have 
important implications for calculations of star formation 
rates based on dividing the stellar mass of a system by an 
age derived from its oldest stars. If the cluster in question 
had spent a large fraction of its time forming stars very 
slowly because of the feedback effects observed in our 
winds runs, the inferred star formation rate would not be 
representative of the earliest stage of the cluster's evolution. 
More simulations are obviously required to determine the 
likelihood of this being true. 

We note that the morphology of the gas in our 
wind-influenced simulations prevents the formation of 
more-or-less spherical wind bubbles. In our simulations, 
the high relative velocities of the wind sources to each other 
and to the gas, and the continuing large scale accretion 
flows present in the model clusters prevented the formation 
of such structures. It is actually very difficult to tell which 
of our simulations include winds and which do not from 
the gas morphology alone, except that the calculations 
including winds exhibit larger quantities of dense shocked 
gas. It may be true therefore that it is not always easy to 
see the effects of winds on young (i.e < lMyr old) clusters. 

In this work, we have not included the effects of pho- 
toionising radiation. The wind sources in our calculations 
are massive enough to have strong ionising fluxes and 
these would also influence the evolution of the protocluster. 
It is not clear whether winds and photoionisation acting 
together help or hinder each other in unbinding stellar 
systems - winds may help clear out dense material from 
around the massive stars, allowing ionising radiation to 
escape, but swept up shells of wind material may confine 
the HII regions at relatively small radii and prevent them 
having much effect. We leave this interesting question to a 
later paper. 
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Figure 16. Column density maps projected along the z-axis of the run without feedback (left panel), the run with isotropic winds 
(centre panel) and the run with anisotropic winds (right panel). White dots represent stars with the large dots denoting stars whose 
masses exceed 10 Mq (i.e. the stars which are wind sources, or those that would be if winds were included in the run without feedback). 
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Figure 17. Evolution with time of the ratio of the mass expulsion rate to the star formation rate (both measured in units of solar masses 
per initial global freefall time) in the control run (left panel), in the isotropic winds run (centre panel), and in the anisotropic winds 
run(left panel). 



